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' Abstract. The time dependent Ginzburg-Landau equation for a two-dimensional granular shear flow is numerically solved, 

_ where we study both the transient dynamics and the steady state of the order parameter. The structural changes of the numerical 

solutions are qualitatively similar to the shear bands observed in the discrete element method (DEM) simulation of the two- 
dimensional granular shear flow. 
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INTRODUCTION 

Flows of granular particles have been well studied due to the importance in technology, engineering, geophysics, 
astrophysics, applied mathematics and physics [1, 2, 3, 4]. The characteristic properties of granular flows are mainly 
caused by the inelastic collisions between particles [5]. Among various studies of granular flows, the study of the 
granular gases under a plane shear plays an important role from many aspects, e.g., the application of the kinetic 
theory to granular gases [6, 7, 8, 9, 10, 11, 12, 13, 14, 15], the shear band formation in moderately dense granular 
gases [16, 17], the long-time tail and the long-range correlation [18, 19, 20, 21, 22, 23, 24, 25, 26, 27], the pattern 
formation in dense granular flow [28, 29, 30, 31, 32, 33], the determination of the constitutive equation for dense 
granular flow [34, 35, 36], as well as the jamming transition [37, 38, 39, 40, 41, 42, 43]. 

The granular hydrodynamic equations derived by the kinetic theory well describe the dynamics of moderately dense 
granular gases [8, 9, 10, 11, 12, 13, 14, 15], though the validity of the kinetic theory is questionable in the case of 
granular gases, because of the lack of scale separation and the existence of the long range correlations, etc [4]. For a 
granular shear flow, a homogeneous state is unstable in the presence of the plane shear [44, 45, 46, 47, 48, 49, 50]. As 
a result, two shear bands are formed near the boundary, and they collide to form one shear band in the center region 
under a physical boundary condition [16, 17]. A similar shear band formation is also observed under the Lees-Edwards 
boundary condition. It is known that both the transient dynamics and the steady state of the hydrodynamic fields can 
be approximately reproduced by the granular hydrodynamic equations [17]. 

To understand the shear band formation after the homogeneous state becomes unstable, we have to develop a 
weakly nonlinear analysis. Recently, Shukla and Alam carried out a weakly nonlinear analysis of granular shear 
flow, where they derived the Stuart-Landau equation of the order parameter defined as the amplitude of disturbance to 
the hydrodynamic fields under a physical boundary condition starting from a set of granular hydrodynamic equations 
[51, 52, 53]. They found the existence of subcritical bifurcations in both relatively dilute and dense regions, while 
a supercritical bifurcation appears in the moderate density region. The Stuart-Landau equation, however, does not 
include any spatial degrees of freedom and cannot be used to study the time evolution of the shear band. 

It is also notable that Khain found the coexistence of a solid phase and a liquid phase in the molecular dynamics 
simulation of a dense granular shear flow [32, 33]. He also demonstrated the hysteresis of the order parameter defined 
as the difference of the densities between the boundary and the center region. It should be noted, however, the 
mechanism of the subcritical bifurcation based on a set of hydrodynamic equations differs from that observed in 
the jamming transition of frictional particles [43]. 

In our previous work, we have developed the weakly nonlinear analysis of a two-dimensional granular shear flow 
and derived the time dependent Ginzburg-Landau (TDGL) equation of the order parameter defined as the amplitude 
of disturbance to the hydrodynamic fields under the Lees-Edwards boundary condition [54]. We introduced a hybrid 
approach to the weakly nonlinear analysis and the resultant TDGL equation is a two-dimensional partial differential 
equation associated with the time dependent diffusion coefficients [54]. The TDGL equation derived by the hybrid 



approach is useful to understand the structural changes of shear bands and we also discussed the bifurcation of the 
order parameter. However, we have not analyzed the solution of the TDGL equation and compared the solution with 
the DEM simulation yet. 

In this paper, we numerically solve the TDGL equation derived in Ref. [54] to exhibit the transient dynamics and 
the steady state of the order parameter. In the following, we review our previous work of the weakly nonlinear analysis 
at first. At second, we present the numerical solution of the TDGL equation. Finally, we discuss and conclude our 
results. 



OVERVIEW OF WEAKLY NONLINEAR ANALYSIS 

In this section, we review our previous results of the weakly nonlinear analysis [54]. At first, we introduce the 
hydrodynamic equations of the area fraction, the velocity fields and the granular temperature. At second, we derive the 
one-dimensional TDGL equation by the ordinary weakly nonlinear analysis. At third, we derive the two-dimensional 
TDGL equation by adopting the hybrid approach to the weakly nonlinear analysis. 



Basic Equations 

Let us introduce our setup and basic equations. We adopt the Lees-Edwards boundary condition for the boundary 
of a two-dimensional granular shear flow, where the upper and the lower image cells move to the opposite directions 
with a constant speed U /2, and the distance between the upper and the lower image cells is given by L [55]. Since 
we assume the two-dimensional granular disks are identical, the mass, the diameter and the restitution coefficient of 
granular disks are given by m, d and e, respectively. In the following argument, we scale the mass, the length and the 
time by m, d and 2d/U, respectively. Therefore, the shear rate U /Lis reduced to e = 2d/L in our units and e is a small 
parameter in the hydrodynamic limit L^> d. 

We employ a set of granular hydrodynamic equations derived by Jenkins and Richman [14]. Although their original 
equations include the angular momentum and the spin temperature, the spin effects are localized near the boundary 
[56] and the effect of rotation can be absorbed in the normal restitution coefficient, if the friction constant is small 
[57, 58]. Therefore, we neglect the rotational degrees of freedom and the dimensionless hydrodynamic equations are 
given by 



(d, +v-V)v = -vV-v 
v(<9,+v-V)v = -VP 
(v/2)(<3 ; +v-V)0 = -P:Vv 
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where v, v = (u,w), 6, t and V = (d/d x ,d/d y ) are the area fraction, the dimensionless velocity fields, the dimension- 
less granular temperature, the dimensionless time and the dimensionless gradient, respectively. The pressure tensor 
P = (P/j), the heat flux q and the energy dissipation rate % are given by 



Pij = 
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respectively, where p(v)6, t,{v)9 l l 2 , r}(v)9 l l 2 , k(v)6^ 2 and A(v)0 3 / 2 are the dimensionless forms of the static 
pressure, the bulk viscosity, the shear viscosity, the heat conductivity and the coefficient associated with the density 
gradient, respectively, and ey = (V^v; + VjVj — 5, y V • v)/2 = x,y) is the deviatoric part of the strain rate. The 
explicit forms of them are listed in Table 1 , where 



*(v) = 



l-7v/T6 
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is the radial distribution function at contact which is only valid for v < 0.7 [59, 60, 61, 62]. 



TABLE 1. The functions in Eqs.(4)-(6). 
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A set of homogeneous solutions of Eqs. (l)-(3) is readily found as 0o = (vo, £y,0, Oq), where Vo and Oq « e 2 / (1 — e 2 ) 
are the mean area fraction and the mean granular temperature, respectively. In our analysis, 9q ~ 0(1). Therefore, 
£ ~ y/l — e 2 and the small e corresponds to the small inelasticity, where e is close to unity [54]. 



Weakly Nonlinear Analysis 

The homogeneous solution 0o is linearly unstable and the disturbance to the hydrodynamic fields with the most 
unstable mode develops as time goes on [44, 45, 46, 47, 48, 49, 50]. To understand the time evolution of 0, we need to 
carry out a weakly nonlinear analysis. For this purpose, we introduce the long time scale T = E 2 t and the long length 
scales (§, Q = e(x,y), respectively. Then, the neutral solution is given by 

n =A L (C,T)0L e ^+c.c. , (8) 

where c.c. represents the complex conjugate and 0^ is the Fourier coefficient of the most unstable mode q c = (0,q c ). 
The amplitude A l (£,t) is independent on because any modes in the sheared frame q(r) = (q^,q^ — £tq^) with 
q% ^ are linearly stable [44, 45, 46, 47, 48, 49, 50]. 
We expand A (£, z) into the series of £ as 

A L (C, T) = £AV(C, t) + £ 2 A!r(C, T ) + £ 3 A^(C, T) + . . . (9) 

and substitute Eqs. (8) and (9) into the hydrodynamic equations (l)-(3). Collecting each order terms of £, we find the 
first non-trivial equation of A\(£, t) at 0(e i ), which is the TDGL equation 

d T A\ = o c A\ +Dd\A\ + PA\\A\\ 2 , (10) 

where D and j3 are the functions of Vo (listed in Table 2 of Ref. [54]) and a c is the maximum growth rate at q c scaled 
by £ 2 . Because of the scaling relations D = D and j3 = £j3, Eq. (10) is rewritten as the equation of the scaled amplitude 

A\{t;,x) = eV*A\{S,x) 

d T A\ = o c A\+bdlA\ + pA\\A\\ 2 . (11) 

The solution of Eq. (11) converges only if ^3 < 0, i.e., in the case of a supercritical bifurcation. 
Developing a similar procedure till 0(e 5 ), we obtain the higher order equation 

d r A L = a c A L + DdjA L + pA L \A L \ 2 + efA L \A L \ 4 + 0(e 3 ) , (12) 

where A L (£, t) = e^ 2 [A\(^, t) + £A^(£, t) + e 2 A\{^, t)] and y is the function of v (listed in Table 2 of Ref. [54]). 
If 7 < 0, the solution of Eq. (12) converges even if j3 > 0, i.e., in the case of a subcritical bifurcation. 



Hybrid Approach to the Weakly Nonlinear Analysis 

The amplitudes A\(£, t) and A L (£, t) are independent of t, and cannot describe the two-dimensional structure of 
shear bands. Thus, we need to introduce a new approach to the weakly nonlinear analysis to derive the two-dimensional 
TDGL equation and study the shear band formation in the granular shear flow. 



At first, we add a small deviation to the most unstable mode as q(f) = q c + Sq(z) and assume n is unchanged if 
the deviation Sq(z) is small 

^ n ^A L (^C,T)^ e ">W- z + c.c, (13) 

where we introduced z = and the amplitude A l (£,£,t) also depends on Combining the contribution from 

the linearly stable mode ^ with n , we introduce the hybrid solution 

Ik = {A L (^C,T)< c +A NL (^C,T)^ ) } e '- ( 'W- z + c.c. 

~ A($,C,T){# e + *g ) } e *W* + c.c., (14) 

where A NL (^ , z) and are the amplitude and the Fourier coefficient, respectively, and we have used a strong 

assumption that A L (^ , z) and A nl (^ , £,T ) are scaled by the common amplitude A , t). Because any modes q(f) 
with q'^ 7^ are linearly stable, decays to zero in the long time limit [44, 45, 46, 47, 48, 49, 50]. 
If we carry out the weakly nonlinear analysis by expanding A(£, , z) as 

A(^C,T)= £ A 1 (^C,T) + e 2 A 2 (^C,T)+ e 3 A 3 (^C^) + ... , (15) 

and using 0h instead of n , we find the two-dimensional TDGL equation 

d z Ai =(7 c Ai+Di{T)dlA l +D 2 (T)d^Ai+Dd^Ai+pAi\Ai\ 2 (16) 

of the rescaled amplitude Ai(|,£,t) = e 1 ' 2 Ai(^,^,z) at (9(e 3 ), where £>i(f) and £>2( T ) are tne ti me dependent 
diffusion coefficients (given by Eqs. (64) and (65) in Ref. [54]). Similarly, we also find the higher order equation 

d r A = a c A+Di(T)dlA+D 2 (T)d i dt ; A+Dd 2 A + pA\A\ 2 + eYA\M 4 + 0(e 3 ) (17) 

Of A(§ , C,T) = £ 1 /2 {Al , Cj T) + eA2 (^ ( £ T ) + £ 2 A3 (^ , C? T)} . 

The two-dimensional TDGL equations (16) and (17) can be solved in the cases of the supercritical bifurcation and 
the subcritical bifurcation, respectively. Because the time dependent diffusion coefficients 5i(t) and Di{x) are the 
functions of ty^-y they decay to zero as time goes on. Therefore, Eqs. (16) and (17) are respectively reduced to Eqs. 
(11) and (12) in the long time limit. 



NUMERICAL SOLUTION OF THE TDGL EQUATION 

In this section, we numerically solve the two-dimensional TDGL equations, where we find the transient dynamics and 
the steady state of the solutions are qualitatively similar to the evolution of the shear band in the two-dimensional 
granular shear flows. 



Numerical Method 

To solve Eqs. (16) and (17) numerically, we prepare the L* x L* square box with the dimensionless system size 
L* = L/d and divide the system into the 100 x 100 grids. We discretize Ai(^,£,t) and A(§,£,t) as A;j(Tfc) and 
A/ j(ffc), respectively, where the continuous variables are given by £, = ix d%, £ = j x dC, and T = k x dz with the small 
increments dE, =dC, = L*/100 and dz = 1.0 x 10~ 4 . We adopt the fourth order Runge-Kutta method to integrate the 
time derivatives and the central difference method to calculate the diffusion terms, e.g., the diffusion terms of A, j(Zk) 
are discretized as 

,2i / _ \ A i+lJ (z k )-2Aij(z k ) +A,_ l .j(z k ) 

"£ A ij{ T k) = ' ' °> 

d^A.M = ^+i( T *)-^^-'(^-^-^( T *)+^-^-'( T *), ( i9) 
d 2 Aij{Tk) = AwW-TAM+Atj-te) ^ (20) 








FIGURE 1. Solutions of Eqs. (16) and (17), where (a-1), (a-2) and (a-3) show the time evolution of Ai (£, £, t) for the supercritical 
bifurcation (v () = 0.20), and (b-1), (b-2) and (b-3) show the time evolution of A(£ , £, t) for the subcritical bifurcation (v () = 0.26), 
respectively. Here, (a-1) and (b-1) are the results of T = 0.70, (a-2) and (b-2) are the results of t = 1.15, and (a-3) and (b-3) are the 
results of T = 1.35, respectively. 



respectively. Since we adopt the Lees-Edwards boundary condition in the weakly nonlinear analysis, we solve Eqs. 
(16) and (17) under the periodic boundary conditions in the sheared frame. The initial values A; ,-(0) and Au(0) are 
given by the superpositions of the sine functions sin(iKW^ + jK^dQ with the random wave numbers Kt and K^, 



Results 

In our weakly nonlinear analysis, the coefficients j3 and y are determined by the mean area fraction Vo [54]. If 
Vo < 0.245, j3 < and the solution of Eq. (16), i.e., A\{£,, £, t), converges. In this case, the supercritical bifurcation of 
the steady amplitude is expected. If 0.245 < V < 0.275, /3 >0 and y< 0. Thus, the solution of Eq. (17), i.e., A(<fj , £, t), 
converges and the subcritical bifurcation of the steady amplitude is expected. Unfortunately,^ > and y > in the 
dense regime Vo > 0.275 and neither Eqs. (16) nor (17) can be used. In the following, we use the small parameter 
e = 0.01 and show the numerical solutions of Eqs. (16) and (17) with Vo = 0.20 and 0.26, respectively. 

Figure 1 displays the numerical solutions of the two-dimensional TDGL equations, where (a-1), (a-2) and (a-3) 
are the time evolution of Ai(£,£,t), and (b-1), (b-2) and (b-3) are the time evolution of A(4,£,t), respectively. In 
both cases, the disturbance in the short wave length is suppressed in the early stage (Figs, l(a-l) and (b-1)) and the 
disturbance in the long wave length survives (Figs. l(a-2) and (b-2)). Then, the shear band is generated in the center 
of the system (Figs. l(a-3) and (b-3)). The steady amplitudes are homogeneous in the ^-direction, because the time 
dependent diffusion coefficients D\ (t) and £>2(t) disappear in the long time limit [54]. As can be seen, we cannot find 
any significant differences between Ai , £, t) and A(<fj , £, t). 

Figure 2 displays the numerical solutions averaged over the ^-direction in the supercritical (Fig. 2(a)) and the 
subcritical (Fig. 2(b)) regimes, respectively. Figure 3 displays the steady amplitudes, where the shear bands have 
peaks at the center of the system £ = 0. 




c c 

FIGURE 2. Time evolution of the order parameter as a function of £ for (a) supercritical (Vo = 0.20) and (b) subcritical 
(Vo = 0.26) regimes, respectively. The open squares, the open circles and the open triangles are the results of T = 0.70, 1.15 
and 1.35, respectively. 
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FIGURE 3. Steady amplitudes as functions of f . The open circles and the closed circles represent |Ai| for Vo = 0.20 and |A| for 
Vq = 0.26, respectively. 



These results are qualitatively similar to the previous result of the area fraction obtained by the DEM simulation [17]. 
The detailed comparison between the DEM simulation and our analysis presented here will be reported elsewhere. 



DISCUSSION AND CONCLUSION 

We numerically solved the two-dimensional TDGL equations (16) and (17) obtained by the hybrid approach to the 
weakly nonlinear analysis, where the structural evolution and the steady state of the solutions are qualitatively similar 
to the shear band observed in the DEM simulation. We also confirmed that the disturbance in the short wave length 
is suppressed in the early stage, and the shear band is survived in the longest wave length. The steady amplitudes are 
homogeneous in the sheared direction, which corresponds to the absence of the time dependent diffusion coefficients 
D\ (t) and £>2(t) in the long time limit. Neither the one-dimensional TDGL equation nor the Stuart-Landau equation 
[51, 52, 53] cannot reproduce such a structural evolution of shear band. 

In conclusion, the solutions of the two-dimensional TDGL equations reproduce the evolution of shear band, which 
are similar to that observed in the DEM simulation [17]. 
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